<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<html>
<head>
<meta http-equiv="Content-Type" content="text/html;charset=iso-8859-1">
<title>MRtrix #VERSION# documentation</title>
<link rel="stylesheet" href="../stylesheet.css" type="text/css" media=screen>
</head>
<body>

<table class=nav>
  <tr>
    <td><a href="dwi.html"><img src="../left.png"></a></td>
    <td><a href="index.html"><img src="../up.png"></a></td>
    <td><a href="../index.html"><img src="../home.png"></a></td>
    <th>Data pre-processing</th>
    <td><a href="roi.html"><img src="../right.png"></a></td>
  </tr>
</table>

<p>
Here, we describe the basic pre-processing steps required for tractography using constrained spherical deconvolution.
We will also show how to generate a number of useful additional images and maps, whose use is recommended.
</p>

<h2><a name='mask'>Brain mask</a></h2>
<p>
A brain mask image is useful for a number of reasons. 
First, it can be used to clean up tensor maps by removing noisy voxels outside the brain.
Second, it can be used to prevent tracks generated by the fibre-tracking algorithm from propagating outside the head.
Last, and most important, it can significantly speed up the <a href='#csd'>CSD</a> computation,
by preventing the analysis of non-brain voxels.
</p>
<p>
For the purposes outlined above, the brain mask should have the same dimensions as the DWI data set.
The simplest approach is to extract the <i>b</i>=0 image from the DWI data set (this example assumes that the first volume in the data set is a <i>b</i>=0 volume):
</p>
<pre>
&gt; <b><a href='../commands/average.html'>average</a> dwi.mif -axis 3 - | <a href='../commands/threshold.html'>threshold</a> - - | <a href='../commands/median3D.html'>median3D</a> - - | <a href='../commands/median3D.html'>median3D</a> - mask.mif</b>
<a href='../commands/average.html'>average</a>: averaging... 100%
<a href='../commands/threshold.html'>threshold</a>: finding min/max... 100%
<a href='../commands/threshold.html'>threshold</a>: building histogram... 100%
<a href='../commands/threshold.html'>threshold</a>: thresholding at intensity 25.5115... 100%
<a href='../commands/median3D.html'>median3D</a>: median filtering... 100%
<a href='../commands/median3D.html'>median3D</a>: median filtering... 100%
</pre>
<p>
<strong>Note:</strong> this example uses <em>pipes</em> to execute a chain of commands.
These are described in more detail <a href='../general/cmdline.html#pipes'>here</a>.
</p>
<p>
The threshold value is in this case determined using simple histogram binning.
You can specify your own threshold if you wish, by supplying the appropriate option to <kbd>threshold</kbd>.
</p>
<p> While this procedure should produce a reasonable mask of the brain, it is
by no means fail-proof. <b>You are strongly encouraged to check the brain mask
using <a href='../general/mrview.html'>MRView</a></b>.  If you need to edit the
mask (for example, if some parts of the brain were not included), you can do so
using the <a href='../general/mrview.html#roi'>ROI analysis</a> sidebar tool
within <a href='../general/mrview.html'>MRView</a>.
</p>

<p class=sep><a href="#top">top</a></p>
<h2><a name='tensor'>Diffusion tensor images</a></h2>
<p>
This section describes the steps necessary to the diffusion tensor and associated parameters 
(see e.g. <a href='../appendix/refs.html#basser'>Basser & Jones, 2002</a> for a review).
The various tensor maps are generated as follows:
</p>
<h3>Tensor components:</h3>
<pre>
&gt; <b><a href='../commands/dwi2tensor.html'>dwi2tensor</a> dwi.mif dt.mif</b>
<a href='../commands/dwi2tensor.html'>dwi2tensor</a>: converting DW images to tensor image... 100%
</pre>
<p>
<strong>Note:</strong> in case the DWI data set headers do not include the DW scheme,
you will need to supply your own encoding file (see <a href='dwi.html#dwscheme'>here</a> for details). 
This can be done as follows (assuming <kbd>encoding.b</kbd> is the encoding file):
<pre>
&gt; <b><a href='../commands/dwi2tensor.html'>dwi2tensor</a> dwi.mif -grad encoding.b dt.mif</b>
<a href='../commands/dwi2tensor.html'>dwi2tensor</a>: converting DW images to tensor image... 100%
</pre>

<h3>Fractional anisotropy (FA) map:</h3>
<pre> 
&gt; <b><a href='../commands/tensor2FA.html'>tensor2FA</a> dt.mif - | <a href='../commands/mrmult.html'>mrmult</a> - mask.mif fa.mif</b>
<a href='../commands/tensor2FA.html'>tensor2FA</a>: generating fractional anisotropy map... 100%
<a href='../commands/mrmult.html'>mrmult</a>: multiplying... 100%
</pre>
<p>
Note that the FA map generated above has been multiplied by the brain mask image to remove the noisy background.
</p>

<h3>Eigenvector (EV) map:</h3>
<pre> 
&gt; <b><a href='../commands/tensor2vector.html'>tensor2vector</a> dt.mif - | <a href='../commands/mrmult.html'>mrmult</a> - fa.mif ev.mif</b>
<a href='../commands/tensor2vector.html'>tensor2vector</a>: generating major eigenvector map... 100%
<a href='../commands/mrmult.html'>mrmult</a>: multiplying... 100%
</pre>
<p>
Note that the EV map generated above has been scaled by the FA map.
</p>


<p class=sep><a href="#top">top</a></p>
<h2><a name='csd'>Constrained spherical deconvolution (CSD)</a></h2>
<p>
This section describes the steps necessary to perform constrained spherical deconvolution 
(for a detailed description of the technique, please refer to <a href='../appendix/refs.html#tournier1'>Tournier <i>et al.</i>. 2004</a> and 
<a href='../appendix/refs.html#tournier2'>Tournier <i>et al.</i> 2007</a>). 
</p>

<h3>Mask of single-fibre voxels:</h3>
First, obtain a mask of high anisotropy voxels. 
These are assumed to contain single fibre-voxels, and will be used to estimate the response function:
</p>
<pre>
&gt; <b><a href='../commands/erode.html'>erode</a> mask.mif -npass 3 - | <a href='../commands/mrmult.html'>mrmult</a> fa.mif - - | <a href='../commands/threshold.html'>threshold</a> - -abs 0.7 sf.mif</b>
<a href='../commands/erode.html'>erode</a>: eroding ... 100%
<a href='../commands/erode.html'>erode</a>: eroding ... 100%
<a href='../commands/erode.html'>erode</a>: eroding ... 100%
<a href='../commands/mrmult.html'>mrmult</a>: multiplying... 100%
<a href='../commands/threshold.html'>threshold</a>: thresholding at intensity 0.7... 100%
</pre>
<p>
The initial erosion step ensure that no edge voxels with artefactually high FA are included in the single-fibre mask.
This mask is then applied to the FA map, and the resulting image is thresholded at FA = 0.7. 
<p>
Note that this value is a guide only - feel free to use a different value if this does not produce satisfactory results. 
Ideally, you should now have a mask containing a few hundred voxels, all
located within high FA white matter regions.  <b>It is very important to check
that the single-fibre mask is suitable, as otherwise the response function
produced in the following step may be totally inappropriate, which would
seriously affect the quality of the CSD output</b>.
If needed, you can edit this mask image to remove unwanted voxels 
using the <a href='../general/mrview.html#roi'>ROI analysis</a> sidebar tool within <a href='../general/mrview.html'>MRview</a>.
</p>

<h3><a name='response'>Response function coefficient:</a></h3>
<p>
The response function SH coefficients can now be estimated from the DW signal in the single-fibre voxels:
</p>
<pre>
&gt; <b><a href='../commands/estimate_response.html'>estimate_response</a> dwi.mif sf.mif response.txt</b>
<a href='../commands/estimate_response.html'>estimate_response</a>: estimating response function... 100%
</pre>
<p>
This should produce the text file <kbd>response.txt</kbd>, containing a series of numbers, such as:
<pre>
232.781 -122.139 52.8959 -17.3446 3.53066 
</pre>
<p>
<strong>Note:</strong> use the <kbd>-grad</kbd> option to supply your own DW scheme if none is to be found in the DWI data set headers.
</p>

<p>
Normally, the values in the response file should alternate between positive and negative while decreasing in magnitude. 
The first value should be positive, and have an amplitude similar to (but not necessarily the same as) the <i>b</i>=0 brain white matter signal 
(you can verify this by loading the <kbd>dwi.mif</kbd> image and placing the focus within a white matter region - 
the intensity is displaced in the bottom right of the statusbar). Note that this is only an approximate guideline: 
it is not unusual for the first coefficient to differ from the b=0 signal quite significantly
(on the other hand, if they differ by orders of magnitude, something has definitely gone wrong).
</p>
<p>
You can also use the <kbd><a href='../commands/disp_profile.html'>disp_profile</a></kbd> command to display the response function:
<pre>
&gt; <b><a href='../commands/disp_profile.html'>disp_profile</a> -response response.txt</b>
</pre>
<p>
This will open a window with a 3D rendering of the response function:
</p>
<img src='response_function.png'>
<p>
The response function should be broadest in the axial plane, and have low amplitude along the <i>z</i>-axis.
If this is not the case, it is likely that the <kbd>estimate_response</kbd> program has used a maximum harmonic order 
that is too high given the data. You can override this value using the <kbd>-lmax</kbd> option, for example:
</p>
<pre>
&gt; <b><a href='../commands/estimate_response.html'>estimate_response</a> dwi.mif sf.mif -lmax 6 response.txt</b>
<a href='../commands/estimate_response.html'>estimate_response</a>: estimating response function... 100%
</pre>
<p>
In general, a value of 6 or 8 is adequate (the number provided <strong>must</strong> be even). 
As a guide, this number should be chosen such that 
the number of distinct DW directions used in the acquisition is greater than the number of parameters that need to be estimated.
For convenience, a table of the number of parameters required for various maximum harmonic orders is provided below:
</p>
<table class=centered>
  <tr><th>Maximum harmonic order (lmax)</th><th>Number of parameters required</th></tr>
  <tr><td>2</td><td>6</td></tr>
  <tr><td>4</td><td>15</td></tr>
  <tr><td>6</td><td>28</td></tr>
  <tr><td>8</td><td>45</td></tr>
  <tr><td>10</td><td>66</td></tr>
  <tr><td>12</td><td>91</td></tr>
  <tr><td><i>n</i></td><td>&frac12; (<i>n</i>+1)(<i>n</i>+2)</td></tr>
</table>


<h3>CSD computation:</h3>
<p>
The CSD computation itself can now be performed:
</p>
<pre>
&gt; <b><a href='../commands/csdeconv.html'>csdeconv</a> dwi.mif response.txt -lmax 8 -mask mask.mif CSD8.mif</b>
<a href='../commands/csdeconv.html'>csdeconv</a>: performing constrained spherical deconvolution... 100%
</pre>
<p>
Here, the maximum harmonic order lmax has been set to 8, which should be 
suitable for most DWI data sets.  You may need to alter this value for 
different acquisition parameters.  Also, the brain mask has been used here to 
prevent unnecessary computations in non-brain voxels; this can speed up the 
computation time significantly.
</p>
<p>
On multi-core systems, the computation time can also be reduced significantly using parallel processing.
The <kbd><a href='../commands/csdeconv.html'>csdeconv</a></kbd> command is capable of running in multi-threaded mode.
To use this feature, set the <kbd>NumberOfThreads</kbd> parameter in the <a href='../appendix/config.html'>MRtrix configuration file</a>.
</p>
<strong>Note:</strong> use the <kbd>-grad</kbd> option to supply your own DW scheme if none is to be found in the DWI data set headers.

<p>
Once the processing is done, you can display the results using the <a href='../general/mrview.html#orientation_plot'>orientation plot</a>
sidebar tool within <a href='../general/mrview.html'>MRview</a>.
</p>


<p class=sep><a href="#top">top</a></p>
<table class=nav>
  <tr>
    <td><a href="dwi.html"><img src="../left.png"></a></td>
    <td><a href="index.html"><img src="../up.png"></a></td>
    <td><a href="../index.html"><img src="../home.png"></a></td>
    <th><a href='#top'>top</a></th>
    <td><a href="roi.html"><img src="../right.png"></a></td>
  </tr>
</table>

<p class=footer>
Donald Tournier<br>
MRtrix version #VERSION#<br>
Last updated #MTIME#
</p>

</body>
</html>


